## Deserving Government Assistance? 
## Public Support for Aid to Struggling Firms and Workers
#############################################################
## Chris Witko
## Tobias Heinrich
#############################################################
## 5) Heterogeneous treatment effects
#############################################################


## Load prepped survey data
###########################
data <- read.csv(file="output/Prepped_survey_data.csv")
CCES <- read.csv(file="output/Prepped_CCES_data.csv")

translate_reason <- data.frame(reason=c("Laid off due to domestic competition", "Laid off due to foreign competition",
                                        "Operations were moved to a foreign country", 
                                        "Laid off due to company mismanagement", "Fired for poor performance", 
                                        "Job was replaced by technology (like robot, artificial intelligence)",
                                        "Laid off due to economic recession", "Company was closed due to pandemic public health restrictions",
                                        "Economic recession", "Competition due to superior use of technology (like robots, artificial intelligence)",
                                        "Mismanagement", "Foreign competition due to lower labor cost",
                                        "Operations closed due to pandemic public health restrictions"),
                               prettyreason=c("Domestic competition", "Foreign competition", "Offshoring",  
                                              "Company mismanagement", "Poor performance",
                                              "Technology", "Recession", "Public health",
                                              "Recession", "Technology", "Mismanagement", "Foreign competiton", "Public health"))

data$W_reasonforjobloss2 <- mapvalues(data$W_reasonforjobloss,
                                      from=translate_reason$reason,
                                      to=translate_reason$prettyreason)
data$F_reasonforstruggle2 <- mapvalues(data$F_reasonforstruggle,
                                       from=translate_reason$reason,
                                       to=translate_reason$prettyreason)

CCES$id <- sample(1:10, size=nrow(CCES), replace=TRUE)
eb_form <- id ~ 1 + Gender_female + Age + Education_uni + I(Party == "Democrat") + I(Party == "Republican") +
  + I(Job == "Full-time") + I(Job == "Unemployed")
X_cces <- model.matrix(eb_form, CCES)[, -1]


## Workers: Calculate heterogeneous effects
###########################################
## Workers: Calculate marginal means
u_salary <- c("$20,000", "$160,000")
u_reasonforjobloss <- na.omit(unique(data$W_reasonforjobloss2))
mmeans <- array(NA, dim=c(1000, length(u_salary), length(u_reasonforjobloss)),
                dimnames=list(Iter=NULL, Salary=u_salary, Reason=u_reasonforjobloss))
for(i in 1:length(u_reasonforjobloss))
{
  for(j in 1:length(u_salary))
  {
    this_data <- subset(data, Which_experiment == "Worker" & W_salary == u_salary[j] &
                          W_reasonforjobloss2 == u_reasonforjobloss[i])
    
    X_data <- model.matrix(eb_form, this_data)[, -1]
    X_eb <- rbind(X_data, X_cces)
    eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
    eb <- ebalance.trim(eb)
    w <- eb$w / sum(eb$w) * nrow(this_data)
    mod <- lm(formula=Choice ~ 1, data=this_data, weights=w)
    VCV <- vcovBS(x=mod, cluster=this_data$id)
    m <- coef(mod)[1]
    se <- sqrt(VCV)[1,1]
    mmeans[, j, i] <- rnorm(1000, m, se)
  }
}


## Workers: Reason | Salary
all_het_worker_reason <- c()
for(i in 1:dim(mmeans)[3])
{
  for(j in 1:dim(mmeans)[3])
  {
    tmp <- adply(.data=mmeans[, , i] - mmeans[, , j], .margins=2,
                 .fun=function(x) data.frame(Mean=mean(x), Q025=quantile(x, 0.025),
                                             Q975=quantile(x, 0.975)))
    tmp$From <- u_reasonforjobloss[j]
    tmp$To <- u_reasonforjobloss[i]
    tmp$tag <- paste0(sort(c(i, j)), collapse="")
    all_het_worker_reason <- rbind(all_het_worker_reason, tmp)
  }
}
all_het_worker_reason <- subset(all_het_worker_reason, Mean != 0)
all_het_worker_reason$Change <- paste0(all_het_worker_reason$From, " --> ", all_het_worker_reason$To)

all_het_worker_reason <- ddply(.data=all_het_worker_reason, .variables=c("tag"),
                               .fun=function(x) subset(x, Change == x$Change[which.max(x$Mean)]))

tmp <- subset(all_het_worker_reason, Salary == "$20,000")
all_het_worker_reason$Change <- factor(all_het_worker_reason$Change, levels=tmp$Change[order(tmp$Mean)])
g <- ggplot(data=all_het_worker_reason, aes(x=Change, y=Mean, ymin=Q025, ymax=Q975,
                                            group=Salary, colour=Salary))
g <- g + geom_hline(yintercept=0, linewidth=0.3)
g <- g + geom_point(size=3, position=position_dodge(width=0.2))
g <- g + geom_linerange(linewidth=0.8, position=position_dodge(width=0.2))
g <- g + coord_flip() + xlab("") + ylab("AMCE")
g <- g + scale_colour_manual(values=c("black", "gray70")) + theme_minimal()
g <- g + theme(axis.text.y=element_text(hjust=0, size=rel(1.34)),
               axis.text.x=element_text(size=rel(1.34)),
               legend.position=c(0.9, 0.14))
ggsave(file="output/Het_worker_reason_salary.pdf", width=10, heigh=6, device=cairo_pdf)


## Workers: Salary | Reason
all_het_worker_reason <- adply(.data=1:8, .margins=1, 
                               .fun=function(x) data.frame(Reason=u_reasonforjobloss[x], Mean=mean(mmeans[, 2, x] - mmeans[, 1, x]), 
                                                           Q025=quantile(mmeans[, 2, x] - mmeans[, 1, x], 0.025),
                                                           Q975=quantile(mmeans[, 2, x] - mmeans[, 1, x], 0.975)))
all_het_worker_reason$Reason <- factor(all_het_worker_reason$Reason,
                                       levels=all_het_worker_reason$Reason[order(all_het_worker_reason$Mean)])
g <- ggplot(data=all_het_worker_reason, aes(x=Reason, y=Mean, ymin=Q025, ymax=Q975))
g <- g + geom_hline(yintercept=0, linewidth=0.3)
g <- g + geom_point(size=3)
g <- g + geom_linerange(linewidth=0.8)
g <- g + coord_flip() + xlab("") + ylab("AMCE")
g <- g + theme_minimal()
g <- g + theme(axis.text.y=element_text(hjust=0, size=rel(1.34)),
               axis.text.x=element_text(size=rel(1.34)))
ggsave(file="output/Het_worker_salary_reason.pdf", width=10, heigh=4, device=cairo_pdf)


## Firms: Calculate heterogeneous effects
#########################################
## Firms: Calculate marginal means: CEO Pay + Salary
u_salary <- c("$20,000", "$70,000", "$160,000")
u_ceopays <- c("$200,000", "$2 million", "$10 million")
mmeans <- array(NA, dim=c(1000, length(u_salary), length(u_ceopays)),
                dimnames=list(Iter=NULL, Salary=u_salary,
                              CEOPay=u_ceopays))
for(i in 1:length(u_ceopays))
{
  for(j in 1:length(u_salary))
  {
    this_data <- subset(data, Which_experiment == "Firm" & F_typicalworkersalary == u_salary[j] &
                          F_ceopay == u_ceopays[i])
    
    X_data <- model.matrix(eb_form, this_data)[, -1]
    X_eb <- rbind(X_data, X_cces)
    eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
    eb <- ebalance.trim(eb)
    w <- eb$w / sum(eb$w) * nrow(this_data)
    mod <- lm(formula=Choice ~ 1, data=this_data, weights=w)
    VCV <- vcovBS(x=mod, cluster=this_data$id)
    m <- coef(mod)[1]
    se <- sqrt(VCV)[1,1]
    mmeans[, j, i] <- rnorm(1000, m, se)
  }
}

all_het_firm_salary_ceo <- adply(.data=1:3, .margins=1, 
                                 .fun=function(x) data.frame(C=u_ceopays[x], Mean=mean(mmeans[, 3, x] - mmeans[, 1, x]), 
                                                             Q025=quantile(mmeans[, 3, x] - mmeans[, 1, x], 0.025),
                                                             Q975=quantile(mmeans[, 3, x] - mmeans[, 1, x], 0.975)))
all_het_firm_ceo_salary <- adply(.data=1:3, .margins=1, 
                                 .fun=function(x) data.frame(C=u_salary[x], Mean=mean(mmeans[, x, 3] - mmeans[, x, 1]), 
                                                             Q025=quantile(mmeans[, x, 3] - mmeans[, x, 1], 0.025),
                                                             Q975=quantile(mmeans[, x, 3] - mmeans[, x, 1], 0.975)))

all_het_firm_salary_ceo$Type <- "Increased typical worker pay"
all_het_firm_salary_ceo$C <- paste0("CEO Pay:\n", all_het_firm_salary_ceo$C)
all_het_firm_ceo_salary$Type <- "Increased CEO pay"
all_het_firm_ceo_salary$C <- paste0("Typical worker pay:\n", all_het_firm_ceo_salary$C)
tmp <- rbind(all_het_firm_salary_ceo, all_het_firm_ceo_salary)

g <- ggplot(data=tmp, aes(x=C, y=Mean, ymax=Q975, ymin=Q025))
g <- g + geom_hline(yintercept=0, linewidth=0.2)
g <- g + geom_point(size=3)
g <- g + facet_wrap(~ Type, ncol=1, scales="free_x")
g <- g + xlab("Conditional") + ylab("AMCE")
g <- g + geom_linerange(linewidth=0.8) + theme_minimal()
g <- g + theme(axis.text.y=element_text(hjust=0, size=rel(1.34)),
               axis.text.x=element_text(size=rel(1.34)))
ggsave(file="output/Het_firm_salary_ceopay.pdf", width=10, height=4, device=cairo_pdf)


## Firms: Calculate marginal means: Reason + Salary
u_salary <- c("$20,000", "$160,000")
u_reasonforstruggle <- na.omit(unique(data$F_reasonforstruggle2))
mmeans <- array(NA, dim=c(1000, length(u_salary), length(u_reasonforstruggle)),
                dimnames=list(Iter=NULL, Salary=u_salary,
                              Reason=u_reasonforstruggle))
for(i in 1:length(u_reasonforstruggle))
{
  for(j in 1:length(u_salary))
  {
    this_data <- subset(data, Which_experiment == "Firm" & F_typicalworkersalary == u_salary[j] &
                          F_reasonforstruggle2 == u_reasonforstruggle[i])
    
    X_data <- model.matrix(eb_form, this_data)[, -1]
    X_eb <- rbind(X_data, X_cces)
    eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
    eb <- ebalance.trim(eb)
    w <- eb$w / sum(eb$w) * nrow(this_data)
    mod <- lm(formula=Choice ~ 1, data=this_data, weights=w)
    VCV <- vcovBS(x=mod, cluster=this_data$id)
    m <- coef(mod)[1]
    se <- sqrt(VCV)[1,1]
    mmeans[, j, i] <- rnorm(1000, m, se)
  }
}

all_het_firm_salary_reason <- adply(.data=1:length(u_reasonforstruggle), .margins=1, 
                                    .fun=function(x) data.frame(Reason=u_reasonforstruggle[x], Mean=mean(mmeans[, 2, x] - mmeans[, 1, x]), 
                                                                Q025=quantile(mmeans[, 2, x] - mmeans[, 1, x], 0.025),
                                                                Q975=quantile(mmeans[, 2, x] - mmeans[, 1, x], 0.975)))

all_het_firm_reason_salary <- c()
for(i in 1:dim(mmeans)[3])
{
  for(j in 1:dim(mmeans)[3])
  {
    tmp <- adply(.data=mmeans[, , i] - mmeans[, , j], .margins=2,
                 .fun=function(x) data.frame(Mean=mean(x), Q025=quantile(x, 0.025),
                                             Q975=quantile(x, 0.975)))
    tmp$From <- u_reasonforstruggle[j]
    tmp$To <- u_reasonforstruggle[i]
    tmp$tag <- paste0(sort(c(i, j)), collapse="")
    all_het_firm_reason_salary <- rbind(all_het_firm_reason_salary, tmp)
  }
}
all_het_firm_reason_salary <- subset(all_het_firm_reason_salary, Mean != 0)
all_het_firm_reason_salary$Change <- paste0(all_het_firm_reason_salary$From, " --> ", all_het_firm_reason_salary$To)

all_het_firm_reason_salary <- ddply(.data=all_het_firm_reason_salary, .variables=c("tag"),
                                    .fun=function(x) subset(x, Change == x$Change[which.max(x$Mean)]))

tmp <- subset(all_het_firm_reason_salary, Salary == "$20,000")
all_het_firm_reason_salary$Change <- factor(all_het_firm_reason_salary$Change,
                                            levels=tmp$Change[order(tmp$Mean)])
g <- ggplot(data=all_het_firm_reason_salary, aes(x=Change, y=Mean, ymin=Q025, ymax=Q975,
                                                 group=Salary, colour=Salary))
g <- g + geom_hline(yintercept=0, linewidth=0.3)
g <- g + geom_point(size=3, position=position_dodge(width=0.2))
g <- g + geom_linerange(linewidth=0.8, position=position_dodge(width=0.2))
g <- g + coord_flip() + xlab("") + ylab("AMCE")
g <- g + scale_colour_manual(values=c("black", "gray70")) + theme_minimal()
g <- g + theme(legend.position=c(0.9, 0.14))
g <- g + theme(axis.text.y=element_text(hjust=0, size=rel(1.34)),
               axis.text.x=element_text(size=rel(1.34)))
ggsave(file="output/Het_firm_reason_salary.pdf", width=10, height=6, device=cairo_pdf)


## Firms: Calculate marginal means: Reason + CEO Pay
u_ceopays <- c("$200,000", "$10 million")
u_reasonforstruggle <- na.omit(unique(data$F_reasonforstruggle2))
mmeans <- array(NA, dim=c(1000, length(u_salary), length(u_reasonforstruggle)),
                dimnames=list(Iter=NULL, CEOPay=u_ceopays,
                              Reason=u_reasonforstruggle))
for(i in 1:length(u_reasonforstruggle))
{
  for(j in 1:length(u_ceopays))
  {
    this_data <- subset(data, Which_experiment == "Firm" & F_ceopay == u_ceopays[j] &
                          F_reasonforstruggle2 == u_reasonforstruggle[i])
    
    X_data <- model.matrix(eb_form, this_data)[, -1]
    X_eb <- rbind(X_data, X_cces)
    eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
    eb <- ebalance.trim(eb)
    w <- eb$w / sum(eb$w) * nrow(this_data)
    mod <- lm(formula=Choice ~ 1, data=this_data, weights=w)
    VCV <- vcovBS(x=mod, cluster=this_data$id)
    m <- coef(mod)[1]
    se <- sqrt(VCV)[1,1]
    mmeans[, j, i] <- rnorm(1000, m, se)
  }
}

all_het_firm_ceopay_reason <- adply(.data=1:length(u_reasonforstruggle), .margins=1, 
                                    .fun=function(x) data.frame(Reason=u_reasonforstruggle[x], Mean=mean(mmeans[, 2, x] - mmeans[, 1, x]), 
                                                                Q025=quantile(mmeans[, 2, x] - mmeans[, 1, x], 0.025),
                                                                Q975=quantile(mmeans[, 2, x] - mmeans[, 1, x], 0.975)))

all_het_firm_reason_ceopay <- c()
for(i in 1:dim(mmeans)[3])
{
  for(j in 1:dim(mmeans)[3])
  {
    tmp <- adply(.data=mmeans[, , i] - mmeans[, , j], .margins=2,
                 .fun=function(x) data.frame(Mean=mean(x), Q025=quantile(x, 0.025),
                                             Q975=quantile(x, 0.975)))
    tmp$From <- u_reasonforstruggle[j]
    tmp$To <- u_reasonforstruggle[i]
    tmp$tag <- paste0(sort(c(i, j)), collapse="")
    all_het_firm_reason_ceopay <- rbind(all_het_firm_reason_ceopay, tmp)
  }
}
all_het_firm_reason_ceopay <- subset(all_het_firm_reason_ceopay, Mean != 0)
all_het_firm_reason_ceopay$Change <- paste0(all_het_firm_reason_ceopay$From, " --> ", all_het_firm_reason_ceopay$To)

all_het_firm_reason_ceopay <- ddply(.data=all_het_firm_reason_ceopay, .variables=c("tag"),
                                    .fun=function(x) subset(x, Change == x$Change[which.max(x$Mean)]))

tmp <- subset(all_het_firm_reason_ceopay, CEOPay == "$200,000")
all_het_firm_reason_ceopay$Change <- factor(all_het_firm_reason_ceopay$Change,
                                            levels=tmp$Change[order(tmp$Mean)])
g <- ggplot(data=all_het_firm_reason_ceopay, aes(x=Change, y=Mean, ymin=Q025, ymax=Q975,
                                                 group=CEOPay, colour=CEOPay))
g <- g + geom_hline(yintercept=0, linewidth=0.3)
g <- g + geom_point(size=3, position=position_dodge(width=0.2))
g <- g + geom_linerange(linewidth=0.8, position=position_dodge(width=0.2))
g <- g + coord_flip() + xlab("") + ylab("AMCE")
g <- g + scale_colour_manual(values=c("black", "gray70")) + theme_minimal()
g <- g + theme(legend.position=c(0.9, 0.14))
g <- g + theme(axis.text.y=element_text(hjust=0, size=rel(1.34)),
               axis.text.x=element_text(size=rel(1.34)))
ggsave(file="output/Het_firm_reason_ceopay.pdf", width=10, height=6, device=cairo_pdf)


## Garbage collection
#####################
rm(list=ls())


